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We develop a method for investigating nonequilibrium dynamics of an ultracold system that is 
initially at thermal equilibrium. Our procedure is based on the classical fields approximation with 
appropriately prepared initial state. As an application of the method, we investigate the influence of 
thermal fluctuations on the quantum phase transition from an antiferromagnetic to phase separated 
ground state in a spin-1 Bose-Einstein condensate of ultracold atoms. We find that at temperatures 
significantly lower than the critical condensation temperature T„ the scaling law for the number of 
created spin defects remains intact. 
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I. INTRODUCTION 

The ground state phase diagram of an antiferromag¬ 
netic Bose-Einstein condensate was studied experimen¬ 
tally in the regime where spatial and spin degrees of 
freedom are decoupled [l|, demonstrating a phase tran¬ 
sition from an antiferromagnetic phase (where only the 
mp = ±1 Zeeman components are populated) to a 
mixed broken axisymmetric phase (where all three Zee- 
man states can be populated). In this regime the sys¬ 
tem follows the predictions of the single spatial mode ap¬ 
proximation (SMA). In other ranges of parameters, how¬ 
ever, spatial separation of components can occur spon¬ 
taneously, and spin domains may appear in the ground 
state of the system In this case the ground state is 

either the antiferromagnetic phase or is phase separated 
into two domains, and the two possibilities are divided 
by a critical point that is characterized by a critical mag¬ 
netic field. 

A system may become out of equilibrium when it is 
driven through the critical point due to the divergence of 
the relaxation time. If symmetry breaking occurs at the 
same time, this out-of-equilibrium process can produce 
various kinds of defects, depending on the dimensionality 
of the system and the form of the order parameter. The 
unified description of these phenomena is described by 
the Kibble-Zurek mechanism (KZM), which was studied 
in a number of physical systems, from the early Universe 
to ultracold atomic gases [JtIII- Among these, Bose- 
Einstein condensates of ultracold atoms offer realistic 
models of highly controllable and tunable systems [T^ . 
Recent experiment with quasi-one-dimensional ultracold 
atoms confirmed the spontaneous creation of solitons via 
the Kibble-Zurek mechanism (T^ . 

In recent papers 00 , we demonstrated that the 
quantum phase transition from an antiferromagnetic to 
phase separated ground state in a spin-1 Bose-Einstein 
condensate of ultracold atoms exhibits scaling laws 
characteristic for systems displaying universal behav¬ 
ior. Phase separation leads to the formation of spin do¬ 
mains, with the number of domain walls depending on 
the quench time. Interestingly, the Kibble-Zurek scal¬ 


ing law was confirmed only for the dynamics close to the 
critical point. Further evolution led to postselection of 
domains, which gave rise to a second scaling law with a 
different exponent. The postselection was attributed to 
the conservation of an additional quantity, namely the 
condensate magnetization. 

In this paper, we develop a method for investigating 
nonequilibrium dynamics of an ultracold system that is 
initially at thermal equilibrium. We apply this method 
to investigate the effect of nonzero temperature on the 
dynamics of the phase transition and the resulting scal¬ 
ing laws. We consider this problem within the framework 
of the classical theory of a complex field with exact equa- 
tions of motion being the Gross-Pitaevskii equations 0 ]. 
The initial condition for the transition includes both the 
condensate and thermal atoms that introduce thermal 
fluctuations in the system. We sample the initial ther¬ 
mal equilibrium within the Bogoliubov approximation at 
a given temperature and for fixed total number of atoms 
and magnetization. Modes orthogonal to the conden¬ 
sate are thermally populated according to the Bogoli¬ 
ubov transformation. In our work, the field has to be 
interpreted not simply as the condensate wavefunction, 
but rather as the total matter field. We present both the 
results of single realizations of the field, which experimen¬ 
tally correspond to single experimental runs, and results 
of averaging over different initial states. We find that 
while the dynamics of the system can be altered by the 
thermal fluctuations, at relatively low temperatures the 
scaling law for the number of domains in the final state 
is intact. 


II. THE MODEL 

In the following, we consider a dilute, weakly interact¬ 
ing spin-1 Bose-Einstein condensate placed in a homoge¬ 
neous magnetic field pointing along the z axis. For the 
sake of completeness, we recall the details of the model 
used previously in 0, 0 to describe the dynamics of 
a condensate at zero temperature within the truncated 
Wigner approximation. In the following section, we will 
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generalize this approach to the case of a condensate a 
nonzero initial temperature. 

We start with the Hamiltonian H = Hq + Ha, wher 
the symmetric (spin-independent) part is 
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Here the subscripts j = — ,0,-I- denote sublevels wit! 
magnetic quantum numbers along the magnetic field axi 
ruf = —1,0,-|-1, m is the atomic mass, p = Pj - 
is the total atom density, V(x) is the externa 
potential. Here we restricted the model to one dimen 
sion, with the other degrees of freedom confined by 
strong transverse potential with frequency uj±. The spin 
dependent part can be written as 


Ha 
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where Ej are the Zeeman energy levels, the spin density i 
F = where F^^y, z are the spin- 

matrices and Ip = (■!/)_|_, ,'0-) • The spin-independen 
and spin-dependent interaction coefficients are given b 
Co = 2fiuj±{2a2 + ao)/3 > 0 and C2 = 2hujj_{a2 — ao)/3 : 
0, where as is the s-wave scattering length for collidin, 
atoms with total spin S. In the following analytic calcu 
lations we assume the incompressible regime where 


Co » |C 2 | , (3 

which is a good approximation eg. in the case of ®^Rb o 
^^Na spin-1 condensate. 

The total number of atoms N = J pdx and magneti 
zation M = J (p+ — p-) dx are conserved quantities. Ii 
reality, there are processes that can change both N am 
M, but they are relatively weak both in spin-1 ^^Na am 
®^Rb condensates [o, 13 and can be neglected on th 
time scales considered below. 

The linear part of the Zeeman shifts Ej induces a ho 
mogeneous rotation of the spin vector around the di¬ 
rection of the magnetic field. Since the Hamiltonian is 
invariant with respect to such spin rotations, we con¬ 
sider only the effects of the quadratic Zeeman shift Hi- 
For sufficiently weak magnetic field we can approximate 
it by a positive energy shift of the nif = ±1 sublevels 
(5 = (£1+ -I- E- — 2Eq)/2 Ri B^A, where B is the mag¬ 
netic field strength and A = {gi + gj)^/x^/16i?HFS, gi 
and gj are the gyromagnetic ratios of electron and nu¬ 
cleus, /rs is the Bohr magneton, i^nFS is the hyperfine 
energy splitting at zero magnetic field [1,1^. Finally, the 
spin-dependent Hamiltonian Q becomes 


Ha 


dx 


^B‘^{p+ + P-) + y • : 
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FIG. 1. Ground state phase diagram of an antiferromagnetic 
condensate for magnetization M = N/2. The green arrow 
indicates the direction of quench into a phase separated state 
considered in Sec. m 


three phases divided by two critical points at 


Bi — Bq 


M 
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where Bq = C 2 plA and p is the total density. The 

ground state can be (i) antiferromagnetic (2C) with 
Ip = (r/!+,O,'0_) for B < Bi, [ii) phase separated into 
two domains of the 2C and ip = (O,'0o,O) type (po) for 
B € (Bi,B 2 ), or (Hi) phase separated into two domains 
of the pq and ip = (ip+,0,0) type (p+) for B > B 2 
What is more, the antiferromagnetic 2C state remains 
dynamically stable, i.e., it remains a local energy mini¬ 
mum up to a critical field B^ > Bi. Consequently, the 
system driven adiabatically from the 2C phase, across 
the phase boundary Bi , and into the separated phase re¬ 
mains in the initial 2C state up to Be > Bi when the 2C 
state becomes dynamically unstable towards the phase 
separation. 


III. MODELING THERMAL FLUCTUATIONS 
BY THE CLASSICAL FIELDS 

Classical fields method is used for description of ther¬ 
mal effects of a single condensate [3l- Here we extend 
the applicability of the method to the system of spin-1 
Bose-Einstein condensates with antiferromagnetic inter¬ 
actions. 

In the method the condensed and noncondensed parts 
of the j-th component are described by a complex func¬ 
tion ipj. The description takes into account modes k with 
energies lower than the temperature: 

/52 i ,2 

( 6 ) 

2m 

The relation between maximal momentum kmax and the 
temperature is not strictly defined and ricut has no unique 
value 0- In general, results of the classical field are cut¬ 
off dependent. Here, we assume that Ucut = 1/3. 

In our simulations, we consider a lattice model for a 
classical field ipj(x) with lattice spacing dx. We enclose 
the atomic field in the ring-shaped quasi-lD geometry 
with periodic boundary conditions at ±L/2 and V[x) = 
0. The total number of atoms is a constant of motion 


Except for the special cases M = 0,±N, the U(a;) = 0 
ground state phase diagram, shown in Fig. [1] contains 
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as well as the total magnetization 

dx^ (|^/’+(a;)p - |^/>_(a:)p) = M. (8) 

X 


The numerical evolution of the field is governed by dis¬ 
cretized counterpart of the Hamiltonian ([T]), ([2|). 

The initial antiferromagnetic (2C) state at thermal 
equilibrium is prepared by employing Bogoliubov trans¬ 
formation of state 2C: = ii^j -I- Stpj [3, with the 

constrain tjjQ = 0. Linearization of the Gross-Pitaevskii 
equation in small fluctuations Stpj (t, x) around uniform 
background 'ipj decouple fluctuations (5'0± from Sipo- 

Fluctuations 5'ip± are composed of phonon (p) and 
magnon (m) branches 
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with quasiparticle energies 

4^’ = C2PV'Csfc2 [2 (co/c2 ) -k 42P], 
4™^ = C2P\/Csfc2(8n+n_ +i'^k^). 
Normalized modes satisfy 


4'"±.“ = 


(m) , (m) 

< ±n = 




± 1/4 


V2(co/C2) 

N ±1/1 
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( 11 ) 


where n± = p±lp. Here we use the spin healing length 
= h/ kJ2m,C2p. 

The small quadrupole mode fluctuations Sil’a [13 ^i^e 
given by 


5^0= E (12) 


fc=0 


with their gapped spectrum for b < be 


^(0) 


el ' = C 2 P\/[Cs + (1 - b'^W - (1 - bl)^ 


(13) 


and the normalized eigenmodes 

( 52 _& 2 )+^ 2 p 
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Here we use a rescaled dimensionless magnetic field 

B 


b = 


Bo 


(14) 

(15) 


To generate the stochastic initial values of the classi¬ 
cal field we proceed as follows, (i) For each realization, 
we generate the fluctuations 6ijjj at temperature T by in¬ 
troducing thermal population of the Bogoliubov modes. 
In practice we generate complex numbers b^^^ for k ^ 0 
of j = ± components and for all A: of j = 0 component 
according to the probability distribution 

(16) 

TT KB-l 

Here (x) denotes (m), (p) or (0). For a given realiza¬ 
tion, we build up fluctuations Sipj according to Bogoli¬ 
ubov transformations (H]) and m- (ii) Then, we create 
the classical field with the constraint that the total atom 
number and magnetization are fixed. The form of the 
field is il’j = 


a+ 


a_ 


- M - 


(17) 

(18) 


where A^_l = is the number of non condensed 

atoms. 

The ensemble of classical fields created in this way is 
the initial state for the time-dependent Gross-Pitaevskii 
(GP) equations 


.dil^o 

‘nr 


2m 


■ CoP ) '00 


+C 2 [{p+ + P-)0o + 20o'0±0-], 
i_I-± = (- er^n-^AB^ 


dt 


dt 


+C2 [{P+ - P-)'0± + PO0± + '0-'0o] > 
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One of the advantages of the method is the ability to 
calculate various correlation functions in a straightfor¬ 
ward way, taking averages over many realizations of the 
stochastic fields. 

Application of the method is justified for a very low 
temperature only since Bogoliubov transformation is the 
approximate solution. Moreover, very early evolution of 
the GP equations can display some transient effects due 
to the fact that the Bogoliubov transformation used in 
the sampling does not produce an exactly stationary dis¬ 
tribution. We have checked if such transients occur for 
parameters used in our simulations and found them to 
be marginal for the quantities we are interested in. 
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FIG. 2. Momentum distribution of non-condensed fields \Sipj{k)f for mp = 1 (left panel), mp = 0 (central), mp = —1 (right). 
Green: initial distribution; Red: after 100s of equilibration with GPEs, Black: Bogoliubov transformation Eqs. and m- 
Here N = 10^, M = N/2, 6 = 0, with length L — 200 fj.m and uip = 27r x 1000 Hz, other parameters as for sodium. 
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FIG. 3. Spectral density |i/ij(fc, cu)p, for mp = 1 (left panel), mp = 0 (central), mp = —1 (right) after 100s of evolution with 
GPEs. Black dashed lines: quasiparticle energies (unj and m given by the Bogoliubov theory. Parameters of the simulation 
are the same as in fig. [2l 


IV. VALIDITY OF THE BOGOLIUBOV 
TRANSFORMATION 

The dynamics with Gross-Pitaevskii equations leads 
the two component antiferromagnetic initial state to 
equilibrium state for long enough times. The equilibra¬ 
tion given by this coupled set of nonlinear equations is 
surprising because we known they may lead to the co¬ 
herent off-equilibrium spin-mixing dynamics. Indeed, it 
was shown by using mean-field theory and adapting the 
single-spatial-mode approximation, that the condensate 
dynamics is well described by a nonrigid pendulum and 
displays a variety of periodic oscillations [2l|. Fortu¬ 
nately, the period of oscillations depends on the initial 
fraction of atoms in the mp = 0 component and is al¬ 
most infinite in our case. Therefore, the time scale as¬ 
sociated with spin-multimode dynamics is much shorter 
and allows for equilibration of the system. The relaxation 
observed in the system is quite intriguing in particular in 
the context of prethermalization phenomena [^ . 

In the low temperature limit, the Bogoliubov theory 
well describes the equilibrium state of the system. The 
comparison of the result of equilibration with GP equa¬ 
tions to the Bogoliubov theory can be treated as an inde¬ 
pendent test of the last one. To this end, we checked the 
validity of the Bogoliubov transformation, and dn. 


as well as the quasiparticle energies, (fTOl) and (IT^ . We 
start numerical calculations with randomly chosen initial 
fluctuations with flat distribution in momentum space, 
see green points in fig. [21 and superimposed constraints 
of given norm N and magnetization M. Then we let the 
initial state to evolve with Gross-Pitaevskii equations for 
a transient time (100s for results presented in Figs. [2] 
and [3]). Next, we compare the achieved distribution of 
non-condensed field in momentum space dtjjj{k ^ 0) to 
the predictions of the Bogoliubov theory. The compar¬ 
ison shown in Fig. [2] is satisfactory. Analysis in the 
frequency domain allows for a test of the validity of 
the elementary excitation picture, see density plots in 
fig. 121 One can easily recognize the magnon and phonon 
branches for mp = ±1 and the quadrupole branch for 
mp = 0. Once again, the comparison with the Bogoli¬ 
ubov theory (marked by black dashed lines) is satisfac¬ 
tory and shows that the quasiparticle energies m and 
(11311 are recovered by equilibration applied to the initial 
state. 


V. THE KIBBLE-ZUREK MECHANISM 

In [Il[ll| we investigated the phase transition from 
the antiferromagnetic to the phase separated state in a 
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spin-1 Bose-Einstein condensate. This continuous phase 
transition is driven by the change of the magnitude of 
the applied magnetic field. Due to the spatial symme¬ 
try breaking in the phase separated state, the transi¬ 
tion is accompanied by the creation of multiple defects 
in the form of spin domain walls, with the number of do¬ 
mains dependent on the quench time. The concept of the 
Kibble-Zurek mechanism (KZM) relies on the fact that 
the system does not follow the ground state exactly in 
the vicinity of the critical point due to the divergence of 
the relaxation time. The dynamics of the system cease 
to be adiabatic at t ~ —t (here we choose t = 0 in the 
first critical point), when the relaxation time becomes 
comparable to the inverse quench rate 

frei Ri |e/e|, (20) 

where e(t) = B — Be ^ t/rq is the distance of the system 
from the critical point. At this moment, the fluctuations 
approximately freeze, until the relaxation time becomes 
short enough again. After crossing the critical point, dis¬ 
tant parts of the system choose to break the symmetry in 
different ways, which leads to the appearance of multiple 
defects in the form of domain walls between domains of 
2C and po phases. The average number of domains is 
related to the correlation length ^ at the freeze out time 

= ( 21 ) 

where z and iz are the critical exponents determined by 
the scaling of the relaxation time ~ and exci¬ 

tation spectrum uj ~ \k\^, with 2 = 1 in the superfluid. 

Interestingly, the Kibble-Zurek scaling law gives cor¬ 
rect predictions only for the dynamics close to the criti¬ 
cal point. Further on, the post-selection of domains was 
observed, which gave rise to a second scaling law with 
a different exponent. The post-selection was attributed 
to the conservation of an additional quantity, namely the 
condensate magnetization. 

The analytical and numerical calculations of [l3, [13 
were carried out within the zero-temperature limit of the 
truncated Wigner approximation. In this section we use 
the method of Sec. lITTlto estimate the influence of nonzero 
temperature on the dynamics of the phase transition and 
the resulting scaling laws. We find that while the dynam¬ 
ics of the system can be altered by the thermal fluctua¬ 
tions, at relatively low temperatures the scaling law for 
the number of domains in the final state is intact. 

We now describe in detail the scenario of the ex¬ 
periment. The antiferromagnetic spin-1 condensate is 
trapped in a ring-shaped quasi-ID trap with strong trans¬ 
verse confinement and the circumference length L. The 
magnetic field is initially switched off, and the atoms 
are prepared in the homogeneous antiferromagnetic (2C) 
ground state with magnetization set to M = N/2. To 
investigate KZM we increase B linearly as 

Bit) = Bo —, ( 22 ) 



FIG. 4. Formation of spin domains through a modified 
Kibble-Zurek mechanism in a ring-shaped ID geometry with 
ring length L = 200 pm and oj± — 2n x 1000 Hz, for N = 10® 
atoms. The evolution of the density of the m/ = 0 component 
in a single realization of the experiment is shown. The 
top and bottom panels correspond to the zero temperature 
and finite temperature results. The quench time is taken as 
TQ = 100 ms. Total condensed fraction at t = 0 is Nk=o/N = 
0.89. 


to drive the system through the two phase transitions into 
a phase separated state. Then, at t > rq the magnetic 
field is kept constant at the level B = Bo- As described in 
Sec.in the ground state of the system becomes separated 
into the 2C and the po phase at B = Bi, but the initial 
2C state remains metastable until the point B = Be- At 
this point, the system is expected to undergo phase tran¬ 
sition accompanied by spatial symmetry breaking. As 
described above, due to the finite quench time the phase 
transition has a nonequilibrium character, and multiple 
spin domains can develop in the system, instead of two 
as the form of the ground state would suggest. Further, 
at the second critical point B = B 2 , there is no symme¬ 
try breaking accompanying the phase transition and the 
spin-domain landscape remains intact. The mean-field 
critical exponents of the symmetry-breaking phase tran¬ 
sition are z = 1 and v = 1/2, which according to the 
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FIG. 5. Comparison of the averaged number of spin domains 
A^d after the quench as a function of the quench time for 
N — 10® atoms. The points are results of numerical simula¬ 
tions averaged over 100 runs. The square and circle poiirts cor¬ 
respond to zero temperature and finite temperature results. 
The scale is logarithmic on both axes and A'a is decreased by 
two to account for the ground state phase separation into two 
domains. The dashed line is the fit to the power law with 
scaling exponent Ud = —0.73 ± 0.04. The grid size is 2® and 
the cut-off is equal to 1/3. 


formula (|2T]) gives the scaling law for the number of do¬ 
mains as Nd ^ However, as shown in [3, [3|i this 

prediction is correct only for the number domain seeds 
formed close to the critical point. When the domains 
become fully formed, their number is decreased in the 
post-selection process, which is due to the existence of 
an upper limit for the number of domains in a system 
with conserved magnetization M. 

Here, we investigate the influence of the finite tem¬ 
perature on the process described above. In Fig. 0] we 
show the evolution of the density of atoms in the ini¬ 
tially unoccupied m/ = 0 state during the process of 
domain formation. The top figure shows a single real¬ 
ization of the zero-temperature truncated Wigner sim¬ 
ulation, which can be interpreted as a result of a sin¬ 
gle experiment. For comparison, analogous result in the 
case of finite temperature, obtained using the method de¬ 
scribed in Sec. uni is shown in the figure below. While 
this particular example corresponds to a relatively short 
quench time, when the post-selection process is not very 
effective, it is visible that some of the initial fluctuations 
merge to form a single domain instead of two. The effect 
of the finite temperature can be seen as the picture of the 
dynamics being more “fuzzy” in the figure below, but the 
number of created domains and their properties seem to 
be unaffected due to the relatively low temperature. 

This conclusion is further confirmed in Fig.jSJ where we 
show the results of systematic averaging of the number 
of created domains over many realizations of the initial 
distribution. The two sets of points showing the result 
for T = 0 and hnite temperature overlap except for the 


very long times, where some decrease of the number of 
domains in the T > 0 case can be observed. Close in¬ 
spection of the dynamics leads us to account this slight 
decrease (on average by less than one domain) to the 
lower number of domains created at the symmetry phase 
transition when the initial state already contains thermal 
fluctuations. In particular, the finite temperature does 
not affect significantly the Nd ^ Tq scaling law pre¬ 
dicted to result from the domain post-selection process. 
This turns out to be the case for any temperature inves¬ 
tigated by our method based on the Bogoliubov trans¬ 
formation. We note that the scaling law may be affected 
by thermal effects at temperatures higher than the ones 
achievable using the current method. This will be the 
topic of a future study. 

VI. CONCLUSIONS 

In summary, we developed a method for investigating 
nonequilibrium dynamics of an ultracold system that is 
initially at thermal equilibrium. Our procedure is based 
on the classical fields approximation with appropriately 
prepared initial state. We described in details how to 
model thermal fluctuations using the Bogoliubov trans¬ 
formation, and demonstrated its validity by performing 
dynamical equilibration of an antiferromagnetic initial 
state with given total atom number and magnetization. 
We studied the effect of the non-zero temperature on 
the scaling law for the number of domains created via 
the Kibble-Zurek mechanism in antiferromagnetic spinor 
condensates. The effect of the finite temperature is vis¬ 
ible on the evolution of the gas density. We find the 
density to be more ’’fuzzy” but the number of created 
domains and their properties rather unaffected by the 
temperature. This result in a relatively low temperature 
shows the strength of the conservation law, namely the 
conservation of magnetization in our system that enforce 
the scaling law. In addition to the main result, we ob¬ 
served relaxation of the antiferromagnetic initial state to 
a thermal-like equilibrium. This intriguing feature pro¬ 
vides an interesting direction for future work. 
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